Subroutine startup()

  ! Routine contains a lot of necessary but unsightly, largely
  ! preprocessor-based, checks that I'd prefer to keep out of main()
  use ISO_FORTRAN_ENV
  use logs
  use my_mpi
  use myfiles
  use arg_module
  use camera_module
  use particles
  implicit none
  integer*4 :: nargs,iargc,ierr,itype, iarg
  character(len=10) :: CharPE,CharNumPEs

  ! Initialise MPI environment 
  call mpi_init(ierr)
  call mpi_comm_rank(mpi_comm_world, MyPE, ierr)
  call mpi_comm_size(mpi_comm_world, NumPEs, ierr)

  ! Allocate arrays to hold particle numbers per processor
  allocate(NumPart_ByPE(0:NumPEs-1)) 
  allocate(NumGhost_ByPE(0:NumPEs-1))
  NumPart_ByPE(:)  = 0
  NumGhost_ByPE(:) = 0

  ! Open a logfile for each process
#ifdef ONE_LOG
  if(MyPE == 0)then
     logunit = OUTPUT_UNIT
  else
     logunit = MyPE + 100
     open(logunit,file='/dev/null',form='formatted')
  endif
#else
  logunit = MyPE + 100
  write(CharPE,'(I3.3)')MyPE
  logfile = './Log_'//trim(adjustl(CharPE))
  open(logunit,file=trim(logfile),form='formatted')
#endif

  ! Write this to a character for tidyness
  write(CharNumPEs,'(I3.3)')NumPEs

  ! Say hello
  write(logunit,*)
  write(logunit,*)' --------------------------------------------------------------------------'
  write(logunit,*)' GridHDF5 v0.1 :: Rob Crain (ICC, Durham) April 2008'
  write(logunit,*)' Running on '//trim(adjustl(CharNumPEs))//' cores'
  write(logunit,*)' --------------------------------------------------------------------------'
  write(logunit,*)
  call flush(logunit)

  ! Open a logfile for each process
#ifdef ONE_LOG
  write(logunit,*)' -using single log file (stdout)'
#else
  write(logunit,*)' -using one log file per core'
#endif

  ! Report orientation method
#ifdef VELOCITY
#ifdef NO_EOS
#ifdef VEL_IGNORE
  write(logunit,*)' -velocity field will not include EOS particles'
#else
  write(logunit,*)' -velocity field will include EOS particles even though omitted from '
  write(logunit,*)'  smoothed fields'
#endif
#endif

#endif

  ! Report orientation method
#ifdef USE_EIGENVECTORS
  write(logunit,*)' -using eigenvectors to orient simulation'
#else
  write(logunit,*)' -using angles (theta/phi/psi) to orient simulation'
#endif

  ! Check that arguments are correct, read and parse into code variables
#ifndef FORCE_EOS
#ifdef EOS_TEMP
  call abort("code compiled to use EOS_TEMP, but FORCE_EOS has not been enabled")
#endif
#ifdef NO_EOS
  call abort("code compiled to use NO_EOS, but FORCE_EOS has not been enabled")
#endif
#endif

  ! Grab command line arguments
  nargs = iargc()
#ifdef USE_CAM_FILE
#ifdef NONEVOLVING
  ! In this case, we have a static snapshot number 
  ! but an evolving frame number...
  if(nargs /= 3) call abort("usage: main paramfile isnap iframe")
  write(logunit,*)' -acquiring camera params from CameraFile'
  call getarg(1,paramfile)
  call getarg(2,arg_snap)
  call getarg(3,arg_frame)
  read(arg_snap,'(I)')isnap
  read(arg_frame,'(I)')iframe
#else
  ! In this case, snapshot and frame are tied together
  if(nargs /= 2) call abort("usage: main paramfile isnap")
  write(logunit,*)' -acquiring camera params from CameraFile'
  call getarg(1,paramfile)
  call getarg(2,arg_snap)
  read(arg_snap,'(I)')isnap
#endif

! --- end of camerafile ifdef ----
#else
! --- start of non-camerafile ifdef ----

#ifdef USE_EIGENVECTORS
#ifdef VELOCITY
#ifdef NONEVOLVING
  if(nargs /= 19) call abort("usage: main paramfile isnap iframe x y z vx vy vz L xprime yprime zprime")
#else
  if(nargs /= 18) call abort("usage: main paramfile isnap x y z vx vy vz L xprime yprime zprime")
#endif
#else ! VELOCITY
#ifdef NONEVOLVING
  if(nargs /= 16) call abort("usage: main paramfile isnap iframe x y z L xprime yprime zprime")
#else
  if(nargs /= 15) call abort("usage: main paramfile isnap x y z L xprime yprime zprime")
#endif 
#endif ! VELOCITY
#else ! USE_EIGENVECTORS
#ifdef VELOCITY
#ifdef NONEVOLVING
  if(nargs /= 13) call abort("usage: main paramfile isnap iframe x y z vx vy vz L theta phi psi")
#else
  if(nargs /= 12) call abort("usage: main paramfile isnap x y z vx vy vz L theta phi psi")
#endif
#else ! VELOCITY
#ifdef NONEVOLVING
  if(nargs /= 10) call abort("usage: main paramfile isnap iframe x y z L theta phi psi")
#else
  if(nargs /= 9) call abort("usage: main paramfile isnap x y z L theta phi psi")
#endif
#endif ! VELOCITY
! --- end of non-camerafile ifdef ----
#endif 

  write(logunit,*)' -acquiring camera params from command line'
  iarg = 1
  call getarg(iarg,paramfile)    ; iarg = iarg + 1
  call getarg(iarg,arg_snap)     ; iarg = iarg + 1
#ifdef NONEVOLVING               
  call getarg(iarg,arg_frame)    ; iarg = iarg + 1
#endif                           
  call getarg(iarg,arg_x)        ; iarg = iarg + 1
  call getarg(iarg,arg_y)        ; iarg = iarg + 1
  call getarg(iarg,arg_z)        ; iarg = iarg + 1
#ifdef VELOCITY
  call getarg(iarg,arg_vx)        ; iarg = iarg + 1
  call getarg(iarg,arg_vy)        ; iarg = iarg + 1
  call getarg(iarg,arg_vz)        ; iarg = iarg + 1
#endif
  call getarg(iarg,arg_L1)       
  call getarg(iarg,arg_L2)       
  call getarg(iarg,arg_L3)       ; iarg = iarg + 1
#ifdef USE_EIGENVECTORS
  call getarg(iarg,arg_xprime_x) ; iarg = iarg + 1
  call getarg(iarg,arg_xprime_y) ; iarg = iarg + 1
  call getarg(iarg,arg_xprime_z) ; iarg = iarg + 1
  call getarg(iarg,arg_yprime_x) ; iarg = iarg + 1
  call getarg(iarg,arg_yprime_y) ; iarg = iarg + 1
  call getarg(iarg,arg_yprime_z) ; iarg = iarg + 1
  call getarg(iarg,arg_zprime_x) ; iarg = iarg + 1
  call getarg(iarg,arg_zprime_y) ; iarg = iarg + 1
  call getarg(iarg,arg_zprime_z) ; iarg = iarg + 1
#else
  call getarg(iarg,arg_theta)    ; iarg = iarg + 1 
  call getarg(iarg,arg_phi)      ; iarg = iarg + 1
  call getarg(iarg,arg_psi)      ; iarg = iarg + 1
#endif


  read(arg_snap,'(I)')isnap
#ifdef NONEVOLVING               
  read(arg_frame,'(I)')iframe
#endif                           
  read(arg_x,*)x0
  read(arg_y,*)y0
  read(arg_z,*)z0
#ifdef VELOCITY
  read(arg_vx,*)vx0
  read(arg_vy,*)vy0
  read(arg_vz,*)vz0
#endif
  read(arg_L1,*)L(1)
  read(arg_L2,*)L(2)
  read(arg_L3,*)L(3)
#ifdef USE_EIGENVECTORS
  read(arg_xprime_x,*)xprime(1)
  read(arg_xprime_y,*)xprime(2)
  read(arg_xprime_z,*)xprime(3)
  read(arg_yprime_x,*)yprime(1)
  read(arg_yprime_y,*)yprime(2)
  read(arg_yprime_z,*)yprime(3)
  read(arg_zprime_x,*)zprime(1)
  read(arg_zprime_y,*)zprime(2)
  read(arg_zprime_z,*)zprime(3)
#else
  read(arg_theta,*)theta
  read(arg_phi,*)phi
  read(arg_psi,*)psi
#endif
#endif




end subroutine startup
